sd_lst[[3]]<-con08$sd
sd_lst[[4]]<-con09$sd
sd_lst[[5]]<-conboot$sd
ci_lst[[1]]<-con06$ci
ci_lst[[2]]<-con07$ci
ci_lst[[3]]<-con08$ci
ci_lst[[4]]<-con09$ci
ci_lst[[5]]<-conboot$ci
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)[2:s_vec[1]]
time_lst[[2]]<-as.vector(time_07$x)[2:s_vec[2]]
time_lst[[3]]<-as.vector(time_08$x)[2:s_vec[3]]
time_lst[[4]]<-as.vector(time_09$x)[2:s_vec[4]]
time_lst[[5]]<-as.vector(time_boot$x)[2:iter.num]
############
#terminate
sd_sd_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(sd_lst[[j]][1:i])
}
sd_sd_lst[[j]]=temp
}
temp=rep(0,(iter.num-1))
for(i in 2:iter.num){
temp[(i-1)]=sd(sd_lst[[5]][1:i])
}
sd_sd_lst[[5]]=temp
figPath<-"/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/"
pngname=paste(figPath,"SD_lc.pdf",sep="")
my_ylab="SD of Standard Deviation"
plotfun(time_lst,sd_sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,220000),my_ylim=c(0,0.01))
######################################
sd_ci_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(ci_lst[[j]][1:i])
}
sd_ci_lst[[j]]=temp
}
temp=rep(0,iter.num-1)
for(i in 2:iter.num){
temp[(i-1)]=sd(ci_lst[[5]][1:i])
}
sd_ci_lst[[5]]=temp
pngname=paste(figPath,"CI_lc.pdf",sep="")
my_ylab="SD of Absolute CI Width"
plotfun(time_lst,sd_ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,200000),my_ylim=c(0,0.02))
plotfun<-function(x,y,filename,my_ylab,my_xlim,my_ylim){
pdf(file=filename, width = 8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,
cex.lab=2.2,lty=1,lwd=2)
plot(x=x[[1]],y=y[[1]],type="l",
xlab="Time (sec)",ylab=my_ylab,xlim=my_xlim,ylim=my_ylim,col="red")
lines(x=x[[2]],y=y[[2]],type="l",
col="blue")
lines(x=x[[3]],y=y[[3]],type="l",
col="green3")
lines(x=x[[4]],y=y[[4]],type="l",
col="orange")
lines(x=x[[5]],y=y[[5]],type="l",
col="black",lty=2)
legend(max(my_xlim)-8000,max(my_ylim),
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8",
"BLBVS-0.9","BootVS"),
lty=c(1,1,1,1,2),lwd=2,pt.cex=1.5,cex=1.7,
col=c("red","blue","green3","orange","black"))
dev.off()
}
converge<-function(beta_blb,r){
not_zero=apply(beta_blb,2,function(x){sum(x!=0)})/nrow(beta_blb)>0.95
p=ncol(beta_blb)
s=nrow(beta_blb)/r
std.bag=matrix(0,nrow=s,ncol=p)
CI.bag=matrix(0,nrow=s,ncol=p)
for(i in 1:s){
std.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,sd)
CI.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,function(x){quantile(x,0.975)-quantile(x,0.025)})
}
######累计平均标准差
std_cum=apply(std.bag,2,cumsum)
std_avg=std_cum
for(i in 1:nrow(std_cum)){
std_avg[i,]=std_cum[i,]/rep(i,ncol(std_cum))
}
###### 累计置信区间宽度
CI_cum=apply(CI.bag,2,cumsum)
CI_avg=CI_cum
for(i in 1:nrow(std_cum)){
CI_avg[i,]=CI_cum[i,]/rep(i,ncol(CI_cum))
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
converge_boot<-function(beta_boot,r){
not_zero=apply(beta_boot,2,function(x){sum(x!=0)})/nrow(beta_boot)>0.8
p=ncol(beta_boot)
std_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
std_avg[i,j]=sd(beta_boot[1:i,j])
}
}
CI_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
CI_avg[i,j]=quantile(beta_boot[1:i,j],0.975)-quantile(beta_boot[1:i,j],0.025)
}
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
setwd(getwd())
######################################### read data
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
time_06<-read.csv("./time_06.csv")
time_07<-read.csv("./time_07.csv")
time_08<-read.csv("./time_08.csv")
time_09<-read.csv("./time_09.csv")
time_boot<-read.csv("./time_boot.csv")
################################################
################## caculte BLB
r=50
s_vec=c(30,20,10,10)
iter.num=100
#####
sd_lst<-list()
ci_lst<-list()
con06=converge(beta_06,r)
con07=converge(beta_07,r)
con08=converge(beta_08,r)
con09=converge(beta_09,r)
conboot=converge_boot(beta_boot,iter.num)
sd_lst[[1]]<-con06$sd
sd_lst[[2]]<-con07$sd
sd_lst[[3]]<-con08$sd
sd_lst[[4]]<-con09$sd
sd_lst[[5]]<-conboot$sd
ci_lst[[1]]<-con06$ci
ci_lst[[2]]<-con07$ci
ci_lst[[3]]<-con08$ci
ci_lst[[4]]<-con09$ci
ci_lst[[5]]<-conboot$ci
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)[2:s_vec[1]]
time_lst[[2]]<-as.vector(time_07$x)[2:s_vec[2]]
time_lst[[3]]<-as.vector(time_08$x)[2:s_vec[3]]
time_lst[[4]]<-as.vector(time_09$x)[2:s_vec[4]]
time_lst[[5]]<-as.vector(time_boot$x)[2:iter.num]
############
#terminate
sd_sd_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(sd_lst[[j]][1:i])
}
sd_sd_lst[[j]]=temp
}
temp=rep(0,(iter.num-1))
for(i in 2:iter.num){
temp[(i-1)]=sd(sd_lst[[5]][1:i])
}
sd_sd_lst[[5]]=temp
figPath<-"/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/"
pngname=paste(figPath,"SD_lc.pdf",sep="")
my_ylab="SD of Standard Deviation"
plotfun(time_lst,sd_sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,220000),my_ylim=c(0,0.01))
######################################
sd_ci_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(ci_lst[[j]][1:i])
}
sd_ci_lst[[j]]=temp
}
temp=rep(0,iter.num-1)
for(i in 2:iter.num){
temp[(i-1)]=sd(ci_lst[[5]][1:i])
}
sd_ci_lst[[5]]=temp
pngname=paste(figPath,"CI_lc.pdf",sep="")
my_ylab="SD of Absolute CI Width"
plotfun(time_lst,sd_ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,200000),my_ylim=c(0,0.02))
plotfun<-function(x,y,filename,my_ylab,my_xlim,my_ylim){
pdf(file=filename, width = 8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,
cex.lab=2.2,lty=1,lwd=2)
plot(x=x[[1]],y=y[[1]],type="l",
xlab="Time (sec)",ylab=my_ylab,xlim=my_xlim,ylim=my_ylim,col="red")
lines(x=x[[2]],y=y[[2]],type="l",
col="blue")
lines(x=x[[3]],y=y[[3]],type="l",
col="green3")
lines(x=x[[4]],y=y[[4]],type="l",
col="orange")
lines(x=x[[5]],y=y[[5]],type="l",
col="black",lty=2)
legend(max(my_xlim)-80000,max(my_ylim),
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8",
"BLBVS-0.9","BootVS"),
lty=c(1,1,1,1,2),lwd=2,pt.cex=1.5,cex=1.7,
col=c("red","blue","green3","orange","black"))
dev.off()
}
converge<-function(beta_blb,r){
not_zero=apply(beta_blb,2,function(x){sum(x!=0)})/nrow(beta_blb)>0.95
p=ncol(beta_blb)
s=nrow(beta_blb)/r
std.bag=matrix(0,nrow=s,ncol=p)
CI.bag=matrix(0,nrow=s,ncol=p)
for(i in 1:s){
std.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,sd)
CI.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,function(x){quantile(x,0.975)-quantile(x,0.025)})
}
######累计平均标准差
std_cum=apply(std.bag,2,cumsum)
std_avg=std_cum
for(i in 1:nrow(std_cum)){
std_avg[i,]=std_cum[i,]/rep(i,ncol(std_cum))
}
###### 累计置信区间宽度
CI_cum=apply(CI.bag,2,cumsum)
CI_avg=CI_cum
for(i in 1:nrow(std_cum)){
CI_avg[i,]=CI_cum[i,]/rep(i,ncol(CI_cum))
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
converge_boot<-function(beta_boot,r){
not_zero=apply(beta_boot,2,function(x){sum(x!=0)})/nrow(beta_boot)>0.8
p=ncol(beta_boot)
std_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
std_avg[i,j]=sd(beta_boot[1:i,j])
}
}
CI_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
CI_avg[i,j]=quantile(beta_boot[1:i,j],0.975)-quantile(beta_boot[1:i,j],0.025)
}
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
setwd(getwd())
######################################### read data
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
time_06<-read.csv("./time_06.csv")
time_07<-read.csv("./time_07.csv")
time_08<-read.csv("./time_08.csv")
time_09<-read.csv("./time_09.csv")
time_boot<-read.csv("./time_boot.csv")
################################################
################## caculte BLB
r=50
s_vec=c(30,20,10,10)
iter.num=100
#####
sd_lst<-list()
ci_lst<-list()
con06=converge(beta_06,r)
con07=converge(beta_07,r)
con08=converge(beta_08,r)
con09=converge(beta_09,r)
conboot=converge_boot(beta_boot,iter.num)
sd_lst[[1]]<-con06$sd
sd_lst[[2]]<-con07$sd
sd_lst[[3]]<-con08$sd
sd_lst[[4]]<-con09$sd
sd_lst[[5]]<-conboot$sd
ci_lst[[1]]<-con06$ci
ci_lst[[2]]<-con07$ci
ci_lst[[3]]<-con08$ci
ci_lst[[4]]<-con09$ci
ci_lst[[5]]<-conboot$ci
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)[2:s_vec[1]]
time_lst[[2]]<-as.vector(time_07$x)[2:s_vec[2]]
time_lst[[3]]<-as.vector(time_08$x)[2:s_vec[3]]
time_lst[[4]]<-as.vector(time_09$x)[2:s_vec[4]]
time_lst[[5]]<-as.vector(time_boot$x)[2:iter.num]
############
#terminate
sd_sd_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(sd_lst[[j]][1:i])
}
sd_sd_lst[[j]]=temp
}
temp=rep(0,(iter.num-1))
for(i in 2:iter.num){
temp[(i-1)]=sd(sd_lst[[5]][1:i])
}
sd_sd_lst[[5]]=temp
figPath<-"/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/"
pngname=paste(figPath,"SD_lc.pdf",sep="")
my_ylab="SD of Standard Deviation"
plotfun(time_lst,sd_sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,220000),my_ylim=c(0,0.01))
######################################
sd_ci_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(ci_lst[[j]][1:i])
}
sd_ci_lst[[j]]=temp
}
temp=rep(0,iter.num-1)
for(i in 2:iter.num){
temp[(i-1)]=sd(ci_lst[[5]][1:i])
}
sd_ci_lst[[5]]=temp
pngname=paste(figPath,"CI_lc.pdf",sep="")
my_ylab="SD of Absolute CI Width"
plotfun(time_lst,sd_ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,200000),my_ylim=c(0,0.02))
plotfun<-function(x,y,filename,my_ylab,my_xlim,my_ylim){
pdf(file=filename, width = 8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,
cex.lab=2.2,lty=1,lwd=2)
plot(x=x[[1]],y=y[[1]],type="l",
xlab="Time (sec)",ylab=my_ylab,xlim=my_xlim,ylim=my_ylim,col="red")
lines(x=x[[2]],y=y[[2]],type="l",
col="blue")
lines(x=x[[3]],y=y[[3]],type="l",
col="green3")
lines(x=x[[4]],y=y[[4]],type="l",
col="orange")
lines(x=x[[5]],y=y[[5]],type="l",
col="black",lty=2)
legend(max(my_xlim)-90000,max(my_ylim),
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8",
"BLBVS-0.9","BootVS"),
lty=c(1,1,1,1,2),lwd=2,pt.cex=1.5,cex=1.7,
col=c("red","blue","green3","orange","black"))
dev.off()
}
converge<-function(beta_blb,r){
not_zero=apply(beta_blb,2,function(x){sum(x!=0)})/nrow(beta_blb)>0.95
p=ncol(beta_blb)
s=nrow(beta_blb)/r
std.bag=matrix(0,nrow=s,ncol=p)
CI.bag=matrix(0,nrow=s,ncol=p)
for(i in 1:s){
std.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,sd)
CI.bag[i,]=apply(beta_blb[c(1:r)+(i-1)*r,],2,function(x){quantile(x,0.975)-quantile(x,0.025)})
}
######累计平均标准差
std_cum=apply(std.bag,2,cumsum)
std_avg=std_cum
for(i in 1:nrow(std_cum)){
std_avg[i,]=std_cum[i,]/rep(i,ncol(std_cum))
}
###### 累计置信区间宽度
CI_cum=apply(CI.bag,2,cumsum)
CI_avg=CI_cum
for(i in 1:nrow(std_cum)){
CI_avg[i,]=CI_cum[i,]/rep(i,ncol(CI_cum))
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
converge_boot<-function(beta_boot,r){
not_zero=apply(beta_boot,2,function(x){sum(x!=0)})/nrow(beta_boot)>0.8
p=ncol(beta_boot)
std_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
std_avg[i,j]=sd(beta_boot[1:i,j])
}
}
CI_avg=matrix(0,r,p)
for(i in 2:nrow(beta_boot)){
for(j in 1:ncol(beta_boot)){
CI_avg[i,j]=quantile(beta_boot[1:i,j],0.975)-quantile(beta_boot[1:i,j],0.025)
}
}
return(list("sd"=rowMeans(std_avg[,not_zero]),
"ci"=rowMeans(CI_avg[,not_zero])))
}
setwd(getwd())
######################################### read data
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
time_06<-read.csv("./time_06.csv")
time_07<-read.csv("./time_07.csv")
time_08<-read.csv("./time_08.csv")
time_09<-read.csv("./time_09.csv")
time_boot<-read.csv("./time_boot.csv")
################################################
################## caculte BLB
r=50
s_vec=c(30,20,10,10)
iter.num=100
#####
sd_lst<-list()
ci_lst<-list()
con06=converge(beta_06,r)
con07=converge(beta_07,r)
con08=converge(beta_08,r)
con09=converge(beta_09,r)
conboot=converge_boot(beta_boot,iter.num)
sd_lst[[1]]<-con06$sd
sd_lst[[2]]<-con07$sd
sd_lst[[3]]<-con08$sd
sd_lst[[4]]<-con09$sd
sd_lst[[5]]<-conboot$sd
ci_lst[[1]]<-con06$ci
ci_lst[[2]]<-con07$ci
ci_lst[[3]]<-con08$ci
ci_lst[[4]]<-con09$ci
ci_lst[[5]]<-conboot$ci
time_lst<-list()
time_lst[[1]]<-as.vector(time_06$x)[2:s_vec[1]]
time_lst[[2]]<-as.vector(time_07$x)[2:s_vec[2]]
time_lst[[3]]<-as.vector(time_08$x)[2:s_vec[3]]
time_lst[[4]]<-as.vector(time_09$x)[2:s_vec[4]]
time_lst[[5]]<-as.vector(time_boot$x)[2:iter.num]
############
#terminate
sd_sd_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(sd_lst[[j]][1:i])
}
sd_sd_lst[[j]]=temp
}
temp=rep(0,(iter.num-1))
for(i in 2:iter.num){
temp[(i-1)]=sd(sd_lst[[5]][1:i])
}
sd_sd_lst[[5]]=temp
figPath<-"/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/"
pngname=paste(figPath,"SD_lc.pdf",sep="")
my_ylab="SD of Standard Deviation"
plotfun(time_lst,sd_sd_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,220000),my_ylim=c(0,0.01))
######################################
sd_ci_lst=list()
for(j in 1:4){
temp=rep(0,(s_vec[j]-1))
for(i in 2:s_vec[j]){
temp[(i-1)]=sd(ci_lst[[j]][1:i])
}
sd_ci_lst[[j]]=temp
}
temp=rep(0,iter.num-1)
for(i in 2:iter.num){
temp[(i-1)]=sd(ci_lst[[5]][1:i])
}
sd_ci_lst[[5]]=temp
pngname=paste(figPath,"CI_lc.pdf",sep="")
my_ylab="SD of Absolute CI Width"
plotfun(time_lst,sd_ci_lst,pngname,my_ylab=my_ylab,
my_xlim=c(0,200000),my_ylim=c(0,0.02))
setwd(getwd())
beta_06<-read.csv("./beta_06.csv")
beta_07<-read.csv("./beta_07.csv")
beta_08<-read.csv("./beta_08.csv")
beta_09<-read.csv("./beta_09.csv")
beta_boot<-read.csv("./beta_boot.csv")
p<-ncol(beta_06)
prop_06<-sort(apply(beta_06,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_07<-sort(apply(beta_07,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_08<-sort(apply(beta_08,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_09<-sort(apply(beta_09,2,function(x)sum(x!=0)/length(x)),decreasing=T)
prop_boot<-sort(apply(beta_boot,2,function(x)sum(x!=0)/length(x)),decreasing=T)
####选择cut-off
pdf(file="/Users/zhangzhang/Desktop/BLB_reversion/jds_blb_reversion/fig/Cutoff_lc.pdf", width =8,height = 6, bg = "white")
par(mar=c(4,5,2,2),mfrow=c(1,1),cex.axis=2.2,cex.lab=2.2,lty=1,lwd=2)
plot(c(1:p),prop_06,
type="o",ylab="Selection Proportion",xlab="Number of Variables",
col="red",ylim=c(0,1),pch=1)
lines(c(1:p),prop_07,
type="o",col="blue",pch=1)
lines(c(1:p),prop_08,
type="o",col="green",pch=1)
lines(c(1:p),prop_09,
type="o",col="orange",pch=1)
lines(c(1:p),prop_boot,
type="o",col="black",pch=1)
legend(1,0.9,
c("BLBVS-0.6","BLBVS-0.7","BLBVS-0.8","BLBVS-0.9","BootVS"),
col=c("red","blue","green","orange","black"),
lwd=2,pt.cex=1.4,cex=1.2)
dev.off()
